Here we compare AR(3) versus TV-persistence forecasts for 24 hours at h={1,2,6,12,24} hour-ahead forecasts separately
Possible improvements:
using RCall
R"library(tvReg)"
using CSV, DataFrames,GLM
using Distributions, LinearAlgebra, Statistics
using Plots, StatsBase, StatsPlots, Colors
using BSON
using LinearAlgebra
using PrettyTables
using Dates
# load main functions
include("tvPersistence_functions_electricity.jl");
myrainbow=reverse(cgrad(:RdYlBu_7, 7, categorical = true));
I take the hourly data for the first hour
Plot of raw data
data_read=CSV.File("DE1523.csv",missingstring=["NA"],header=true) |> DataFrame;
datesall=Date.(first.(data_read.Column1,10),dateformat"y-m-d");
datesallyr=Date.(first.(data_read.Column1,4),dateformat"y");
v=first.(data_read.Column1,4)
pos=findall(x->x==1,v[2:end].!==v[1:end-1]);
v=first.(data_read.Column1[1:24:end],4)
poshour=findall(x->x==1,v[2:end].!==v[1:end-1]);
plot(data_read.Price,size=(1000,200),title="Hourly prices all",label=false,frame=:box,color=myrainbow[1],titlefontsize=10,xticks=([1;pos],2015:2023))
plot!(fontfamily="serif-roman",titlefontsize=10, xtickfontsize=10,ytickfontsize=10,ylabelfontsize=10)
# transform data
datatransf=[]
a=[]
b=[]
for i=1:24
datatransf0,a0,b0 = transformation_asinh(data_read.Price[i:24:end]);
push!(datatransf,datatransf0)
push!(a,a0)
push!(b,b0)
end
PLots of transformed data we work with, shaded area is in-sample (rolling window OOS is the rest)
for i in 1:24
plot(datatransf[i],size=(1000,200),title="Prices - hour $(i-1)",label=false,frame=:box,color=myrainbow[1],
xticks=([1;poshour],2015:2023),titlefontsize=10)
plot!(fontfamily="serif-roman",titlefontsize=10, xtickfontsize=10,ytickfontsize=10,ylabelfontsize=10)
display(plot!([1:1024],fill=myrainbow[1], fillalpha=0.2,linecolor=:transparent,seriestype=:vspan,label=false))
end
Here we have ratios TV-EWD to AR(3) losses in terms of RMSE and MAE, so that 0.9 is 10% gain of TV-EWD across all 24 hours at all horizons
folder="results_transformed" # transformed data not tuned (not transformed back)
#folder="results" # transformed data not tunted kernel
#folder="Old_results" # nontransformed data
forecasts1=[]
for i in 1:24
forecast_0 = BSON.load("./$folder/results_electricity_$i-1.bson")
push!(forecasts1,forecast_0["forecasts"])
end
forecasts2=[]
for i in 1:24
forecast_0 = BSON.load("./$folder/results_electricity_$i-2.bson")
push!(forecasts2,forecast_0["forecasts"])
end
forecasts6=[]
for i in 1:24
forecast_0 = BSON.load("./$folder/results_electricity_$i-6.bson")
push!(forecasts6,forecast_0["forecasts"])
end
forecasts12=[]
for i in 1:24
forecast_0 = BSON.load("./$folder/results_electricity_$i-12.bson")
push!(forecasts12,forecast_0["forecasts"])
end
forecasts24=[]
for i in 1:24
forecast_0 = BSON.load("./$folder/results_electricity_$i-24.bson")
push!(forecasts24,forecast_0["forecasts"])
end
losses1=[[rmse(forecasts1[i][:,1],forecasts1[i][:,2]) for i=1:24]'
[rmse(forecasts1[i][:,1],forecasts1[i][:,4]) for i=1:24]'
[mae(forecasts1[i][:,1],forecasts1[i][:,2]) for i=1:24]'
[mae(forecasts1[i][:,1],forecasts1[i][:,4]) for i=1:24]'];
losses2=[[rmse(forecasts2[i][:,1],forecasts2[i][:,2]) for i=1:24]'
[rmse(forecasts2[i][:,1],forecasts2[i][:,4]) for i=1:24]'
[mae(forecasts2[i][:,1],forecasts2[i][:,2]) for i=1:24]'
[mae(forecasts2[i][:,1],forecasts2[i][:,4]) for i=1:24]'];
losses6=[[rmse(forecasts6[i][:,1],forecasts6[i][:,2]) for i=1:24]'
[rmse(forecasts6[i][:,1],forecasts6[i][:,4]) for i=1:24]'
[mae(forecasts6[i][:,1],forecasts6[i][:,2]) for i=1:24]'
[mae(forecasts6[i][:,1],forecasts6[i][:,4]) for i=1:24]'];
losses12=[[rmse(forecasts12[i][:,1],forecasts12[i][:,2]) for i=1:24]'
[rmse(forecasts12[i][:,1],forecasts12[i][:,4]) for i=1:24]'
[mae(forecasts12[i][:,1],forecasts12[i][:,2]) for i=1:24]'
[mae(forecasts12[i][:,1],forecasts12[i][:,4]) for i=1:24]'];
losses24=[[rmse(forecasts24[i][:,1],forecasts24[i][:,2]) for i=1:24]'
[rmse(forecasts24[i][:,1],forecasts24[i][:,4]) for i=1:24]'
[mae(forecasts24[i][:,1],forecasts24[i][:,2]) for i=1:24]'
[mae(forecasts24[i][:,1],forecasts24[i][:,4]) for i=1:24]'];
a=plot([losses1[2,:]./losses1[1,:] losses2[2,:]./losses2[1,:] losses6[2,:]./losses6[1,:] losses12[2,:]./losses12[1,:] losses24[2,:]./losses24[1,:]],
ylim=[0.5,1.2],ylab="TV-EWD / AR",xlab="hours",label=["h=1" "h=2" "h=6" "h=12" "h=24"],
frame=:box,title="RMSE Ratio for all hours and horizons",xticks=(1:24,0:23),
color=[myrainbow[1] myrainbow[3] myrainbow[5] myrainbow[6] myrainbow[7]],linewidth=1)
a=hline!([1],label=false,color=:black)
b=plot([losses1[4,:]./losses1[3,:] losses2[4,:]./losses2[3,:] losses6[4,:]./losses6[3,:] losses12[4,:]./losses12[3,:] losses24[4,:]./losses24[3,:]],
ylim=[0.5,1.2],ylab="TV-EWD / AR",xlab="hours",label=["h=1" "h=2" "h=6" "h=12" "h=24"],
frame=:box,title="MAE Ratio for all hours and horizons",xticks=(1:24,0:23),color=[myrainbow[1] myrainbow[3] myrainbow[5] myrainbow[6] myrainbow[7]],
linewidth=1)
b=hline!([1],label=false,color=:black)
plot(a,b,layout=(1,2),size=(1100,400))
plot!(fontfamily="serif-roman",titlefontsize=10, xtickfontsize=10,ytickfontsize=10,ylabelfontsize=10)
println("***************************************************************")
println(" RMSE - h=1 MAE - h=1")
println("***************************************************************")
display([["Hours" "AR" "TV-EWD" "AR" "TV-EWD"];
[[i for i=0:23]'
losses1]'])
println("***************************************************************")
println(" RMSE - h=2 MAE - h=2")
println("***************************************************************")
display([["Hours" "AR" "TV-EWD" "AR" "TV-EWD"];
[[i for i=0:23]'
losses2]'])
println("***************************************************************")
println(" RMSE - h=6 MAE - h=6")
println("***************************************************************")
display([["Hours" "AR" "TV-EWD" "AR" "TV-EWD"];
[[i for i=0:23]'
losses6]'])
println("***************************************************************")
println(" RMSE - h=12 MAE - h=12")
println("***************************************************************")
display([["Hours" "AR" "TV-EWD" "AR" "TV-EWD"];
[[i for i=0:23]'
losses12]'])
println("***************************************************************")
println(" RMSE - h=24 MAE - h=24")
println("***************************************************************")
display([["Hours" "AR" "TV-EWD" "AR" "TV-EWD"];
[[i for i=0:23]'
losses24]'])
pockets of predictability identified for all hours, those in darker blue area are not spuriously (or randomly) generated, we use simple persistent process in the bootstrap
Areas simply point to periods when TV-EWD model improves accuracy over AR(3)
poshouroos=poshour.-1024;
poshouroos=poshouroos[3:end];
bw=0.03;
for i in 1:24
display(plot_pockets(forecasts1[i][:,2].-forecasts1[i][:,1],forecasts1[i][:,4].-forecasts1[i][:,1],bw,"Pockets of predictability - hour $(i-1)"))
end
Lubos's comparison period
naive is 9,2559 LEAR 4.4655
[mean([mae(forecasts1[i][616:616+554-1,1],forecasts1[i][616:616+554-1,2]) for i=1:24])
mean([mae(forecasts1[i][616:616+554-1,1],forecasts1[i][616:616+554-1,4]) for i=1:24])]
v=first.(data_read.Column1[9:24:end][1024:end],4)
poshouroos=findall(x->x==1,v[2:end].!==v[1:end-1]);
plot(forecasts1[9][:,[1;2;4]],size=(1200,220),label=["data" "AR forecast" "TV-EWD forecast"],frame=:box,title="h=1 forecasts",
color=[myrainbow[1] myrainbow[5] myrainbow[7]],titlefontsize=10,xticks=(poshouroos,2018:2023))
plot(forecasts6[9][:,[1;2;4]],size=(1200,220),label=["data" "AR forecast" "TV-EWD forecast"],frame=:box,title="h=6 forecasts",
color=[myrainbow[1] myrainbow[5] myrainbow[7]],titlefontsize=10,xticks=(poshouroos,2018:2023))